The goal here is to create a prediction model for the value of benthic over the Chesepeake bay region.
This prediction model will be graphically portrayed with a:
ggplot2 over an actual mapggplot2 over an actual mapplotly surface plotThe following code will automatically install packages that are not already installed and load them into the library:
if(require('pacman')){
library('pacman')
}else{
install.packages('pacman')
library('pacman')
}
Loading required package: pacman
pacman::p_load(sp, EnvStats, gstat, ggplot2, rmarkdown, reshape2, ggmap, RColorBrewer, parallel, dplyr, plotly)
Note that Benthic.df is included in the EnvStats package and is not a data frame but an sp class, which is very similar to a data frame
#Assignments
benthic <- Benthic.df
lat <- benthic$Latitude
lon <- benthic$Longitude
index <- benthic$Index
head(benthic, 3)
The Benthic Data set is an sp class objects and requires for its co-ordinates to be set:
coordinates(benthic) = ~Longitude + Latitude
In order to create a krige prediction model two things are required:
If a 4th degree polynomial model was used to model the benthic values over space, a variogram can be made from that model thusly:
(The quadratic model is a required parameter of the kriging process, in exercise 12.2 a quadratic model will be used, the spacific mathematics of this though is somewhat beyond me, however the kriging method uses the residuals from that surface model in order to make the predictions.)
x <- lon
y <- lat
benthic.vg <- variogram(object = index ~
x + y +
x^2 + x*y + y^2 +
x^3 + x^2*y + x*y^2 + y^3 +
x^4 + x^3*y + x^2*y^2 + x*y^3 + y^4,
data = benthic)
head(benthic.vg)
np dist gamma
1 4545 0.02410565 0.8506164
2 4149 0.06703996 0.9393502
3 3883 0.11396396 1.1513831
4 4023 0.16138943 1.4599238
5 5963 0.20615519 1.5641015
6 6732 0.24983031 1.8992994
dir.hor dir.ver id
1 0 0 var1
2 0 0 var1
3 0 0 var1
4 0 0 var1
5 0 0 var1
6 0 0 var1
benthic.vg.fit.exp <- fit.variogram(benthic.vg, model=vgm(1,"Exp", 0.5,1))
Thus the variogram model needed for the krige prediction has been ascertained
This will essentially act as the resolution of the overlay, Most monitors are about 140 ppi, hence assuming the plot will be a 9 inch square, about 800 x 800 pixels should look right, but 200 x 200 will be the used due to performance.
xgrid <- seq(min(lon), max(lon), length.out = 200)
ygrid <- seq(min(lat), max(lat), length.out = 200)
This domain base however has to be a spatial object, this can be acheived by defining co-ordinates in the data frame
xy.surface <- expand.grid(lon = xgrid, lat = ygrid)
coordinates(xy.surface) = ~lon + lat
head(xy.surface, 3)
SpatialPoints:
lon lat
1 -77.31570 38.0018
2 -77.30859 38.0018
3 -77.30147 38.0018
Coordinate Reference
System (CRS) arguments: NA
benthic.pred <- krige(formula = Index ~1,
locations = benthic,
newdata = xy.surface,
model= benthic.vg.fit.exp)
[using ordinary kriging]
head(benthic.pred)
coordinates var1.pred
1 (-77.3157, 38.0018) 3.144851
2 (-77.30859, 38.0018) 3.140807
3 (-77.30147, 38.0018) 3.136427
4 (-77.29436, 38.0018) 3.131703
5 (-77.28724, 38.0018) 3.126626
6 (-77.28013, 38.0018) 3.121187
var1.var
1 2.321630
2 2.316530
3 2.311443
4 2.306367
5 2.301299
6 2.296237
Now there is a prediction model, if that prediction model is extrapolated over our domain (xy.surface) the various plots can be created
In order to create a heatmap, it is first necessary to take the sp prediction object and make it into pixel data, otherwise the heatmap won’t look quite right:
class(benthic.pred)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
gridded(benthic.pred) = TRUE
class(benthic.pred)
[1] "SpatialPixelsDataFrame"
attr(,"package")
[1] "sp"
spplot(benthic.pred["var1.pred"],main="Predictions of Benthic Index - ordinary kriging")
First it is necessary to create a data frame from the sp object:
benthic.pred.df <- cbind(benthic.pred@coords, benthic.pred@data)
benthic.pred.df.tidy <- melt(benthic.pred.df)
No id variables; using all as measure variables
head(benthic.pred.df)
ggplot() +
geom_tile(data = benthic.pred.df, aes(lon, lat, fill = var1.pred)) +
labs(fill = "Benthic Index") #
This plot would be easier to read if it had discrete intervals rather than a continuous scale:
benthic.pred.df$bin <- cut(benthic.pred.df$var1.pred,
breaks = as.vector(seq(from = 0, to = 5, length.out = 20)))
ggplot() +
geom_tile(data = benthic.pred.df, aes(lon, lat, fill = bin)) +
labs(fill = "Benthic Index",
title = "Heatmap of Benthic Index Values",
x = "Longitude",
y = "Lattitude")
NA
Now by adjusting the colour scale and opacity, this can be used as an overlay for the map from exercise 10.1:
benthic.pred.df$bin <- cut(benthic.pred.df$var1.pred,
breaks = as.vector(seq(from = 0, to = 5, length.out = 9)))
#Get a map
bbox <- make_bbox(benthic.pred.df$lon, benthic.pred.df$lat, f = 0.01)
map <- get_map(location = bbox, maptype = 'toner')
maptype = "toner" is only available with source = "stamen".
resetting to source = "stamen"...
Map from URL : http://tile.stamen.com/toner/9/146/194.png
Map from URL : http://tile.stamen.com/toner/9/147/194.png
Map from URL : http://tile.stamen.com/toner/9/148/194.png
Map from URL : http://tile.stamen.com/toner/9/146/195.png
Map from URL : http://tile.stamen.com/toner/9/147/195.png
Map from URL : http://tile.stamen.com/toner/9/148/195.png
Map from URL : http://tile.stamen.com/toner/9/146/196.png
Map from URL : http://tile.stamen.com/toner/9/147/196.png
Map from URL : http://tile.stamen.com/toner/9/148/196.png
Map from URL : http://tile.stamen.com/toner/9/146/197.png
Map from URL : http://tile.stamen.com/toner/9/147/197.png
Map from URL : http://tile.stamen.com/toner/9/148/197.png
#Create the colours
cols <- brewer.pal(n = length(levels(benthic.pred.df$bin)), name = "Oranges")
#Create the Plot
ggmap(ggmap = map) +
geom_tile(data = benthic.pred.df, aes(lon, lat, fill = bin, alpha = var1.pred))+
scale_fill_manual(values = cols) +
labs(fill = "Benthic Index",
x = "Longitude",
y = "Lattitude",
title = "Heatmap of Benthic Index") +
guides(alpha=FALSE)
NA
In order to view the surface plot, contours can also be used.
Recall from Week 11/Exercise 10, that contour and surface plots work a little differently to ordinary plots.
A domain grid of equally spaced \(x\) and \(y\) values is required and a matrix of \(z\) values that can be overlayed onto that domain grid in order to provide the surface values required to create the surface plot.
The method from Lecture 10/Wk. 11 involved predicting over the square domain of xy.surface, this won’t work here because instead of the predict function, the krige function was used, hence it is necessary to convert the output from the krige function, into the output that would have been given by the predict function, which is the required \(z\) surface matrix.
This can be acheived with the data.matrix function:
z_surface_matrix <- data.matrix(benthic.pred[1])
Using base R packages, a contour plot can be created with the contour function:
contour(xgrid, ygrid, z_surface_matrix,
xlab="-Longitude (degrees West)", ylab="Latitude (degrees North)",
main = "contour plot of Krige Prediction",
labcex = 0.8)
mtext("Krige using a 4th degree polynomial model", cex = 0.5)
In order to create a contour plot in ggplot 2, it will be necessary to use ‘tidy’ data, which is a table of data such that each column corresponds to a variable:
rownames(z_surface_matrix) <- xgrid
colnames(z_surface_matrix) <- ygrid
benthic.pred.melt <- melt(z_surface_matrix,
varnames = c("Longitude", "Latitude"),
value.name = "Benthic_Index")
head(benthic.pred.melt)
Observe that this is another way of getting the data frame that was required in part 1 (benthic.pred.df), it is included here as another method and for the sake of completion, also it matches the method used in week 10/lecture 11.
Now using ggplot2 a contour plot, can be created:
ggplot(data = benthic.pred.melt, aes(x = Longitude,
y = Latitude,
z = Benthic_Index,
colour = ..level..),
size = 0.5) +
labs(x = "Longitude", y = 'Latitude',
title = 'Benthic Index',
subtitle = "Chesapeake Bay, Maryland",
caption = "Higher index value is better",
col = 'Benthic Index', size = 'Benthic Index') +
geom_contour(lwd = 1) +
theme_classic()
The advantage to using ggplot2, is the map overlay can be taken advantage of:
map <- get_map(location = bbox, maptype = 'watercolor')
maptype = "watercolor" is only available with source = "stamen".
resetting to source = "stamen"...
Map from URL : http://tile.stamen.com/watercolor/9/146/194.jpg
Map from URL : http://tile.stamen.com/watercolor/9/147/194.jpg
Map from URL : http://tile.stamen.com/watercolor/9/148/194.jpg
Map from URL : http://tile.stamen.com/watercolor/9/146/195.jpg
Map from URL : http://tile.stamen.com/watercolor/9/147/195.jpg
Map from URL : http://tile.stamen.com/watercolor/9/148/195.jpg
Map from URL : http://tile.stamen.com/watercolor/9/146/196.jpg
Map from URL : http://tile.stamen.com/watercolor/9/147/196.jpg
Map from URL : http://tile.stamen.com/watercolor/9/148/196.jpg
Map from URL : http://tile.stamen.com/watercolor/9/146/197.jpg
Map from URL : http://tile.stamen.com/watercolor/9/147/197.jpg
Map from URL : http://tile.stamen.com/watercolor/9/148/197.jpg
ggmap(ggmap = map) +
labs(fill = "Benthic Index",
x = "Longitude",
y = "Lattitude",
title = "Heatmap of Benthic Index") +
guides(alpha=FALSE) +
geom_contour(data = benthic.pred.melt, aes(x = Longitude,
y = Latitude,
z = Benthic_Index),
col = 'grey40',
binwidth = 0.15,
size = 0.70) +
scale_fill_manual(values = cols)
This can be combined with the heatmap layer:
ggmap(ggmap = map) +
labs(fill = "Benthic Index",
x = "Longitude",
y = "Lattitude",
title = "Heatmap of Benthic Index") +
guides(alpha=FALSE) +
geom_tile(data = benthic.pred.df, aes(lon, lat, fill = bin, alpha = var1.pred)) +
geom_contour(data = benthic.pred.melt, aes(x = Longitude,
y = Latitude,
z = Benthic_Index),
col = 'indianred',
binwidth = 0.2,
size = 0.6) +
scale_fill_manual(values = cols)
A surface plot allows a three dimensional view of the model
In base package, a surface plot can be made:
persp(xgrid, ygrid, z_surface_matrix,
xlim = c(-77.3, -75.9), #the values domain needs to be adjusted above
ylim = c(38.1, 39.5), #Make sure to set an appropriate domain
# zlim = c(0, 6),
theta = -45, phi = 30, d = 0.5,
xlab="-Longitude (degrees West)",
ylab="Latitude (degrees North)",
zlab="Benthic Index", ticktype = "detailed")
surface extends beyond the box
title(main=paste("Surface Plot of Benthic Index", "Loess Smoothing", sep="\n"))
Plotly has the advantage of being fully interactive and not as difficult to constrain within a plot/box area:
p <- plot_ly(x = ygrid, y = xgrid, z = z_surface_matrix) %>%
add_surface() %>%
layout(
title = "Benthic Index in Chesapeake Bay",
scene = list(
xaxis = list(title = "Longitude",
range = c(38.2, 38.9)), #You shouldn't need to edit the
yaxis = list(title = "Latitude", #The domain here, do it above
range = c(-77, -75.9)),
zaxis = list(title = "Benthic Index")
))
p
The surface plot really reinforces the fact that the bay has the lowest (i.e. best) benthic index value.